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Abstract 

The maximum likelihood estiaator has been used 
to extract stability and control derivatives from 
flight data for aany years. Most of the literature 
on aircraft estiaation concentrates on new develop- 
aents and applications, assuaing familiarity with 
basic estiaation concepts. This paper presents 
done of these basic concepts. The paper briefly 
discusses the maximum likelihood estiaator and the 
aircraft equations of motion that the estiaator 
uses. The basic concepts of minimization and esti- 
mation are examined for a simple computed aircraft 
example. Hie cost functions that are to be mini- 
mized during estimation are defined and discussed. 
Graphic representations of the cost functions are 
given to help illustrate the minimization process. 
Finally, the basic concepts are generalized, and 
estimation from flight data is discussed. Some of 
the major cor elusions for the computed example are 
also developed for the analysis of flight data. 

Introduction 

The maximum likelihood estimator has been used 
to obtain stability and control estimates from 
flight data for nearly 20 years. The results of 
many applications have been reported worldwide. 
Reference 1 contains a representative list of some 
of these reports. Several good texts (including 
Refs. 2 and 3) contain thorough treatments of the 
theory of maximum likelihood estimation. Experi- 
ence reports*' 4 ' 5 pointing out practical consid- 
erations for applying the maximum likelihood esti- 
mator have also been published. Stability and 
control derivatives estimated from flight data are 
currently required for correlation studies with 
predictive techniques, handling qualities documen- 
tation, design compliance, aircraft simulator 
enhancement and refinement, and control system 
design. Correlation, simulation, and control sys- 
tem design applications are discussed in Ref. 6. 
Current studies have concentrated on estimation 
model structure determination, 7,5 equation error 
with state reconstruction, 5 '* 0 '* 1 and maximum like- 
lihood estimation in the frequency domain. * 3 

Most of the reports in the estimation area con- 
centrate on new developments and applications, 
assuming familiarity with the basic concepts of 
maximum likelihood estimation. In this paper some 
of these basic concepts are reviewed, concentrating 
on simple, idealized models, These simple models 
provide insights applicable to a wide variety of 
real problems. 

This paper presents some fundamentals of maxi- 
mum likelihood estimation as applied to the air- 
craft problem. It briefly discusses the maximum 
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likelihood estimator and the aircraft equations of 
motion that the estimator uses. The basic aspects 
of minimization and estimation are then examined in 
detail for a simple computed aircraft example. 
Finally, the discussion is expanded to the general 
aircraft estimation problem. 

Symbols 

system matrices 
lateral acceleration, g 

reference span, ft 

coefficient of rolling moment 

coefficient of yawing moment 

coefficient of side-force 

general functions 

measurement noise covariance 
matrix 

acceleration due to gravity, ft/sec 2 

approximation to the information 
matrix 

moment of inertia about subscrip- 
ted axis, slug-ft 2 

l general index 

J cost function 

Kg sidewash factor 

L rolling moment divided by I x , deg/ 

sec 2 , or number of iterations 

l ' rolling moment, ft-lb 

Lyj rolling moment due to yaw jet, 

ft-lb 

m mass, slug 

N number of time points or cases 

n state noise vector or number of 

unknowns 

p roll rate, deg/sec 

q pitch rate, deg/sec 

q dynamic pressure, lb/ft 2 

R innovation covariance matrix 

r yaw rate, deg/sec 


A, B,C, D, F, G 

a y 

b 

Cf 

c n 

Cy 

f (* ) , g( • ) 
GG* 

g 

H 



a 


reference area, ft^ 
t tine, sec 

u control input vector 


{or obtaining maximum likelihood estimates for the 
aircraft problem are given. In the following sec- 
tions, both the concepts and the computations 
involved in a simple but realistic example are dis- 
cussed in detail. 


V forward velocity, ft/sec 

x state vector 

*a <Va > z a distance between lateral accelerom- 

Y V Y eter and the center of gravity 

along the appropriate axis, ft 

z observation vector 

z predicted Kalman-filtered estimate 

a angle of attack, deg 

6 angle of sideslip, deg 

A time sample interval, sec 

6 control deflection, deg 

6 a aileron deflection, deg 

5 e elevon deflection, deg 

6 r rudder deflection, deg 

n measurement noise vector 

0 pitch angle, deg 


The aircraft parameter estimation problem can 
be defined quite simply m general terms, the sys- 
tem investigated is assumed to be modeled by a set 
of dynamic equations containing unknown parameters. 
To determine the values of the unknown parameters, 
the system is excited by a suitable input, and the 
input and actual system response are measured. The 
values of the unknown parameters are then inferred 
based on the requirement that the model response to 
the given input match the actual system response, 
when formulated in this manner, the problem of 
identifying the unknown parameters can be easily 
solved by many methods; however, complicating fac- 
tors arise when application to a real system is 
considered. 

The first complication results from the impos- 
sibility of obtaining perfect measurements of the 
response of any real system. The inevitable sensor 
errors are usually included as additive measurement 
noise in the dynamic model. Once this noise is 
introduced, the theoretical nature of the problem 
changes drastically. It is no longer possible to 
exactly identify the values of the unknown param- 
eters,- instead, the values must be estimated by 
some statistical criterion. The theory of esti- 
mation in the presence of measurement noise is 
relatively straightforward for a system with 
discrete time observations, requiring only basic 
probability. 


U mean 

5 vector of unknowns 

0 standard deviation 

T time, sec 

$ transition matrix or bank angle, 

deg 

fi integral of transition matrix 

Subscripts: 

p,q,r,o,a,0,8, partial derivative with respect 
5,5 a ,6 r ,6 e to subscripted quantity 


The second complication of real systems is the 
presence of state noise. State noise is random 
excitation of the system from unmeasured sources, 
the standard example for the aircraft stability and 
control problem being atmospheric turbulence. If 
state noise is present and measurement noise is 
neglected, the analysis results in the regression 
algorithm. 

When both state and measurement noise are con- 
sidered, the problem is more complex than in the 
cases that have only state noise or only measure- 
ment noise. Reference 14 develops the Maine-Iliff 
formulation of the maximum likelihood estimator in 
continuous/discrete time, which accounts for both 
state and measurement noise. This formulation has 
a continuous system model with discrete sampled 
observations. 


0 b.'as or at time zero 

m measured quantity 

Other nomenclature: 

~ predicted estimate 

estimate 

* transpose 

Maximum Likelihood Estimation 

The concept of maximum likelihood is discussed 
in this section. First the general heuristic prob- 
lem is discussed, and then the specific equations 


The final problem for real systems is modeling. 
It has been assumed throughout the above discussion 
that for some value (called the "correct" value) of 
the unknown parameter vector, the system is cor- 
rectly described by the dynamic model. Physical 
systems are seldom described exactly by simple 
dynamic models, so the question of modeling error 
arises. No comprehensive theorv of modeling error 
is available. The most common approach is to 
ignore it: Any modeling error is simply treated as 

state noise or measurement noise, or both, in spite 
of the fact that the modeling error may be deter- 
ministic rather than random. The assumed noise 
statistics can then be adjusted to include the con- 
tribution of the modeling error. This procedure 
is not rigorously justifiable, but combined with a 
carefully chosen model, it is probably the best 
approach available. 


2 





V 


-i 


i 


With the above diecueeion in Bind, it ie pos- 
eible to Bake a aore preciae, aatheaaticaXly proba- 
bilistic stateaent of the paraaeter estiaation 
problea. the first step is to define the general 
systea aodel (aircraft eguations of aotion). this 
aodsl can be written in the continuous/discrete 
fora as 


X(tQ> - x 0 (1) 

i(t) - f [x(t),u(t),£l + F(£)n(t) (2) 


where 


♦ - exp (A(t i+1 - t 4 > ) 


* 



exp (At) dt B 


When state noise is important, the nonlinear 
fora of Eqs. (1) to (3) is intractable. For the 
linear aodel defined by Bqs. (S) to (7), the cost 
function that accounts for state noise is 


t(tj) - g(x(t i ),u(t i ),C) + Gttlhi (3) 

where x is the state vector, z is the observa- 
tion vector, f and g are systea state and obser- 
vation functions, u is the known control input 
vector, £ is the unknown paraaeter /actor, n ie 
the state noise vector, and n is the aeasureaent 
noise vector. The state noise vector is assumed 
to be zero-mean white Gaussian and stationary, and 
the measurement noise vector is assumed to be a 
sequence of independent Gaussian random variables 
with zero mean and identity covariance. For each 
possible estirvte of the unknown parameters, a 
probability tht the aircraft response tiae his- 
tories attain values near the observed values can 
then be defined. The maximum likelihood estimates 
are defined as those that maximize this probabil- 
ity. Maximum likelihood estimation has many desir- 
able statistical characteristics; for example, it 
yields asymptotically unbiased, consistent, and 

efficient estimates. 1 ^ 

If there is no state noise and the matrix G is 
known, then the maximum likelihood estimator mini- 
mizes the cost function 


J(£> -“-X l*< fc i> - *r<ti>) , R' 1 l*(t i ) 

2 i-1 

- ZjUj^] ♦ ~ N In |R| (11) 

where R is the innovation covariance matrix. The 
zj(t^) term in Eq . (11) is the Kalman-filtered esti- 
mate of z, which, if the state noise covariance is 
zero, reduces to the form of Bq. (4). If there is 
no state noise, the second term of Bq. ( 11 ) is of 
no consequence, (unless one wishes to include 
elements of the G matrix) and R can be replaced by 
GG* which makes Bq. (11) the same as Eq. (4). 

To minimize the cost function J(£), we can 
apply the Newton-Raphson algorithm which chooses 
successive estimates of the vector of unknown 

coefficients, £. Let L be the iteration number. 

The L + 1 estimate of £ is then obtained from the 
L estimate as follows: 

t W 1 - *L - (12) 


J(£) X (Zlti) - Z 5 (t i )]*(GG«)- , (z(t i ) 

2 i-1 

- z^tti)] + ^ N In | (GG*) | (4) 

where GG* is the measurement noise covariance 
matrix, and z^ft^ is the computed response esti- 
mate of z at t^ for a given value of the unknown 
parameter vector £. The cost function is a func- 
tion of the difference between the measured and 
computed time histories. 

If Bqs. (2) and (3) are linearized (as is the 
case for the stability and control derivatives in 


the aircraft problea), 

x(t 0 ) - x 0 (5) 

x(t) - Ax( t) + Bu(t) + Fn(t) (6) 

Z < t ± ) - CxU^ + Du(t^) + GDi (7) 

For the no-state-noise case, the z^(t^) term of 
Eq. (4) can be approximated by 

x^(t 0 ) - x 0 (£) (8) 

X^(t; + i ) - ♦X£(t i ) + ♦(U(t i ) + u(t i+ ,))/2 (9) 

z 5 (ti) - Cx^tt^) + Du( t| ) (10) 


If R is assumed fixed the first and second gra- 
dients are defined as 

N 

V C J(£) = - X l* (t i> - Vv 1 * 

i*1 

(GG*)- 1 (V e z e (ti) 1 (13) 

2 £ ~ , ~ 

V|J(£) - X ! v £ z t (t i )] * (GG * r1lv £ z £ (t i )1 
i=1 

- jhz<ti) - ^<t t )i* 

2 

(GG*)* 1 [V t z e ( t ± ) ) (14a) 

The Gauss-Newton approximation to the second gra- 
dient is 

N 

v|j( C) s X 17 £ z £ (t i ,, * (GG *> _1 [’£Z£( t i )) (14b) 

The Gauss-Newton approximation, which is sometimes 
referred to as modified Newton-Raphson, is com- 
putationally much easier than the Newton-Raphson 
approximation because the second gradient of the 
innovation never needs to be calculated, in addi- 
tion, it can have the advantage of speeding the 
convergence of the algorithm, as is discussed in 
the Simple Aircraft Example section. 
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Figure 1 illustrates the maximum likelihood 
estimation concept. The measured response of the 
aircraft is compared with the estimated response, 
and the difference between these responses is 
called the response error. The cost functions of 
Eqs. ( 4 ) and (11) include this response error. The 
Gauss-Newton computational algorithm is used to 
find the coefficient values that maximize the cost 
function. Each iteration of this algorithm pro- 
vides a new estimate of the unknown coefficients 
on the basis of the response error. These new 
estimates of the coefficients are then used to 
update the mathematical model of the aircraft, pro- 
viding a new estimated response and, therefore, a 
new response error. The updating of the mathemat- 
ical model continues iteratively until a conver- 
gence criterion is satisfied. The estimates result- 
ing from this procedure are the maximum likelihood 
estimates. 

The maximum likelihood estimator also provides 
a measure of the reliability of each estimate based 
on the information obtained from each dynamic 
maneuver. This measure of the reliability, analo- 
gous to the standard deviation, is called the 
Cramer-Rao bound 16 or the uncertainty level. The 
Cramer-Rao bound as computed by current programs 
should generally be used as a measure of relative 
accuracy rather than absolute accuracy. The bound 
is obtained from the approximation of the infor- 
mation matrix, H. This matrix equals the approxi- 
mation to the second gradient given by Bq. (14b). 

The bound for each unknown is proportional to the 
the square root of the corresponding diagonal 
element of H. That is for the 1 th unknown, the 
Cramer-Rao bound is /H(i,i) . 


where 


c ygS + C 


pb „ rb 
+ CV — + 


*p 2V * 1 *r W 


Cy _ 6 + C- 


*0 


c t - c* e e ♦ c tp & ♦ c tr H ♦ c t6 6 


3b 

+ C *0 + C *a It 


c n * C ng® + C n p + c n r 2 v + 


6b 

+ C + c • ■— 
'-no + '-ng 2V 


where the 6 term is summed over all controls. 
The observation equations are 


8m ‘ K e( B - ^ P + —■ r) 


Pm * P 




qs „ 

a * C Y p + — ‘ r 

rm mg 1 g g 


(19) 


( 20 ) 


( 21 ) 


( 22 ) 

(23) 

(24) 

(25) 


y a„ 


(p 2 + r 2 ) 


(26) 


The formulation and minimization algorithm 
discussed above is implemented with the Iliff-Haine 
code (MMLE3 maximum likelihood estimation program). 
The program and computational algorithms are 
described fully in Ref. 17. All the computations 
shown and described in the remainder of this paper 
use the algorithms exactly as described in Ref. 17. 

Aircraft Equations of Motion 

For the discussion that follows in later sec- 
tions of this paper, some knowledge of the aircraft 
equations of motion is assumed. To clarify some of 
that discussion, the aircraft equations are dis- 
cussed briefly in this section. 

Generalized nonlinear equations of motion are 
given in detail in Ref. 17, which fully describes 
the Iliff-Haine code (MMLE3 program). All compu- 
tations and aircraft examples in this paper use 
the linearized form for the lateral-directional 
equations. These equations are given below and 
referred to in the remainder of the paper. 

8 » ^ (C y + 8 0 ) *■ ^ cos 8 sin ♦ 


+ p sin a - r cos a (15) 

pl x " * r xz * 3 sbc i + ^ r(I y - I z> + P^xz (16) 

rlz - P^z = ^ sbc n + P<l(Ix - V * <5 rI xz <>7) 

♦ » p + r cos $ tan 8 

+ q sin $ tan 6 + (18) 


Pm = P + Po ( 27) 

r m » ^ t r 0 (28) 

The state, control, and observation vectors for 
the lateral -directional mode can then be defined as 

x = (8 p r ♦>* (29) 

u « ( 6 a 6 r )» (30) 


z = (6 m Pm r m ♦m a y m Pm r m 1 * 131 * 
Simple Aircraft Example 

The basic concepts involved in a parameter 
estimation problem can be illustrated by using a 
simple example representative of a realistic air- 
craft problem. The example chosen here is repre- 
sentative of an aircraft that exhibits pure roll- 
ing motion from an aileron input. This example, 
although simplified, typifies the motion exhibited 
by many aircraft in particular flight regimes, such 
as the F-14 aircraft flying at high dynamic pres- 
sure, the F— 111 aircraft at moderate speeds with 
the wing in the forward position, and the T-37 
aircraft at low speed. 

Derivation of an equation describing this 
motion is straightforward. Figure 2 shows a sketch 
of an aircraft with the x-axis perpendicular to the 
plane of the figure (positive forward on the air- 
craft). The rolling moment (l"), roll rate (p), 
and aileron deflection (6 a ) are positive as shown. 
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For this example, the only state is p end the only 
control is 6 a . The result of summing moments is 

I x p « L^p, 6 a ) (32) 

The fi'st-order Taylor expansion then becomes 


where 


4 - exp (Lpi) 


4 



exp (LpT ) dr 


H 


L$H - exp (Lpd)] 


p - LpP + L« a « a (33) 

where 

l'- i x l 

Since the aileron is the only control, it is nota- 
tionally simpler to use 6 instead of 6 a for the 
discussion of this example. Equation (33) can then 
be written as 


and A is the length of the sample interval (t ^ +1 
- t^ > • Simplifying the notation 

«i+ 1/2 " < 6 i + 5 i + 1>/ 2 <38> 

then 

Pi +1 ” 4Pi + V^i+i/2 (39) 


p - LpP + L 56 (34) 

An alternate approach that results in the same 
equation is to combine Eq, (16) with Eq. (20). sub- 
stituting for C and then eliminate the terms that 

are zero for our example. This yields 

P*x “ 5 sb ( c t p + C * 6 S ) < 35) 

where p is the roll rate and 6 is the aileron 
deflection. Rearranging terms, the equation can 
be put into the dimensional derivative form of 
Eq. (34). 

Equation (34) is a simple aircraft equation 
where the forcing function is provided by the aile- 
ron and the damping by the damping-in-roll term, L p . 

In subsequent sections we examine in detail the 
parameter estimation problem where Eq. (34) des- 
cribes the system. For this single-degree-of- 
freedom problem, the maximum likelihood estimator 
is used to estimate either Lp or L^ or both for a 

given computed time history. 

We will assume that the system has measurement 
noise, but no state noise as in £qs. ( 1 ), ( 2 ), and 
(3). Equation (4) then gives the cost function 
for maximum likelihood estimation. The weighting 
GG* is unimportant for this problem, so let it 
equal 1. For our example, Eqs. (2) and (3) become 
x A » and Zi = Therefore, Eq. (4) becomes 

1 n 

J<V L «> = j I [ Pi " Pi<V L « )l2 (36) 

i-1 

where Pi is the value of the measured response p at 
time ti and p^(L p ,L{) is the computed time history 
of p at time ti for L p « L p and L 5 = L$. Through- 
out the rest of the paper, where computed data (not 
flight) are used, the measured time history refers 
to pi> and the computed time history refers to 
Pi<Lp,Lj). The computed time history is a function 
of the current estimates of Lp and L$, but the 
measured time history is not. 

The most straightforward method of obtaining 
Pi is with Eqs. ( 8 ) and (9). In terms of the nota- 
tion stated above, 

Pi+i " 4Pi + 4< 6 i + 6 i+i>/2 !37) 


The maximum likelihood estimate is obtained by 
minimizing Eq* (36). The Gauss-Newton method 
described earlier is used for this minimization. 
Equation (12) is used to determine successive 
values of the estimates of the unknowns during the 
minimization. 

For this simp! a problem, £ « [Lp Lj] * and suc- 
cessive estimates of Lp and Lj are determined by 
updating Eq. (12). The first and second gradients 
of Eq. (12) are defined by Eqs. (13) and (14b). 

The complete set of equations are given in Ref. 17. 

The entire procedure can now be written for 
obtaining the maximum likelihood estimates for this 
simple example. To start the algorithm, an initial 
estimate of L p and Lj is needed. This is the value 

AAA 

of Cg. With Eq. (12), Si and subsequently are 

defined by using the first and second gradients of 
J(Lp,Lj) from Eq. (36). The gradients for this 

particular example from Eqs. (13) and (14b) are 

V£o( 5 l ) - - ^ (Pi - Pi>V£ Pi (40) 

i =1 

N 

V|j(S L ) s I (’ 5 Pi>* (V 5 Pi) (41) 

i = 1 

With the specific equations defined in this 
section for this simple example, we ca.i now proceed 
in the next section to the computational details of 
a specific example. 

Computational Details of Minimization 

In the previous section we specified the equa- 
tions for a simple example and described the proce- 
dure for obtaining estimates of the unknowns from a 
dynamic maneuver. In this section we give the com- 
putational details for obtaining the estimates. 

Some of the basic concepts of parameter estimation 
are best shown with computed data where the correct 
answers are known. Therefore, in this section we 
study two examples involving computed time his- 
tories. The first example is based on data that 
have no measurement noise, which results in esti- 
mates that are the same as the correct value. The 
second example contains significant measurement 
noise; consequently, the estimates are not the same 
as the correct values. Throughout the rest of the 
paper, where comouted data is used, the term 
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"no-noise case” is used for the case with no noise 
added and "noisy case" for the case where noise has 
been added. 

Since we are studying a simple computed 
example, it is desirable to keep it simple enough 
to complete some or all of the calculations on a 
home computer or, with some labor, on a calculator. 
With this in mind, the number of data points needs 
to be kept small. For this computed example, 10 
points (time samples) are used. The simulated 
data, which we refer to as the measured data, are 
based on Eq. (34). We use the same correct values 
of L p and Lj (-0.2500 and 10.0, respectively) for 

both examples. In addition, the same input ( 6 ) is 
used for both examples, the sample interval (A) is 
0.2 sec, and the initial conditions are zero. 

Tables of all the significant intermediate values 
are given for each example. These values are given 
to four significant digits, although to obtain 
exactly the same values with a computer or calcula- 
tor requires the use of 13 significant digits, as 
in the computation of these tables. If the four- 
digit numbers are used in the computation, the 
answers will be a few tenths of a percent off, but 
will still serve to illustrate the minimization 
accuracy. In both examples, the initial values of 

Lp and Lg (or ) are -0.5 and 15.0, respectively. 
Example With No Measurement Noise 

The measurement time history for no measurement 
noise (no-noise case) is shown in Fig. 3. The aile- 
ron input starts at zero, goes to a fixed value, 
and then returns to zero. The resulting roll-rate 
time history is also shown. The values of the 
measured roll rate to 13 significant digits are 
given in Table 1 along with the aileron input. 

Table 2 shows the values for Lp, L<j, and J for 
each iteration, along with the values of 4 and 4 
needed for calculations of p^. In three iterations 
the algorithm converges to the correct values to 

four significant digits for both Lp and L<j. L^ 
overshoots slightly on the first iteration and then 

comes quickly to the correct answer. Lp overshoots 
slightly on the second iteration. 

Figure 4 shows the match between the measured 
data and the computed data for each of the first 
three iterations. The match is very good after two 
iterations. The match is nearly exact after three 
iterations. 

Although the algorithm has converged to four- 
digit accuracy in L p and L$, the value of the cost 

function, J, continues to decrease rapidly between 
iterations 3 and 4. This is a consequence of using 
the maximum likelihood estimator on data with no 
measurement noise. Theoretically, using infinite 
accuracy the value of J at the minimum should be 
zero. However, with finite accuracy the value of J 
becomes small but never quite zero. This value is 
a function of the number of significant digits that 
are being used. For the 13-digit accuracy used 
here, the cost eventually decreases to approxi- 
mately 0.3 x 10 ' 28 . 


Example with Measurement Noise 

The data used in this example (noisy case) are 
the same as those used in the previous section, 
except that pseudo-Gaussian noise has been added 
to the roll rate. The time history is shown in 
Fig. 5. Tne signal-to-noise ratio is quite low in 
this example, as is readily apparent by comparing 
Figs. 3 and 5. The exact values of the time his- 
tory to 13-digit accuracy are shown in Table 3. 

The values of Lp, Lj, 4, 4 , and J are shown for 

each iteration in Table 4. The algorithm con- 
verges in four iterations. The behavior of the 
coefficients as they approach convergence is much 
like the no-noise case. The most notable result of 

this case is the converged values of Lp and Lj ( 

which are somewhat different from the correct 
values. The match between the measured and com- 
puted time history is shown in Fig. 6 for each 
iteration. No change in the match is apparent for 
the last two iterations. The match is very good 
considering the amount of measurement noise - 

In Fig. 7, the computed time history for the 
no noise estimates of Lp and Lj is compared to 
that for the noisy-case estimates of Lp and Lj. 

Because the algorithm converged to values somewhat 
different than the correct values, the two com- 
puted time histories are similar but not identical. 

The accuracy of the converged elements can be 
assessed by looking at the Cramer-Rao inequal- 
ity* 8,17 discussed earlier. The Cramer-Rao bound 
can be obtained from the information matrix cor- 
rected for observed noise amplitude as follows. 

11 = £ ^minimum* J* V(N-1) 

The Cramer-Rao bounds for Lp and Lj are the square 
roots of the diagonal elements of the H matrix, or 
/H (1,1) and /H(2,2), respectively. The Cramer-Rao 
bounds are 0.1593 and 1.116 for Lp and Lj, respec- 
tively. The errors in Lp and Lj are less than the 
bounds . 

Cost Functions 

In the previous section we obtained the maximum 
likelihood estimates for computed time histories by * 

minimizing the values of the cost function. To 
fully understand what occurs in this minimization, 
we must study in more detail the form of the cost 
functions and some of their more important charac- 4 

teristics. In this section, the cost function for 
the no-noise case is discussed briefly. The cost 
function of the noisy case is then discussed in 
more detail. The same two time histories studied 
in the previous section are examined here. The 
noisy case is more interesting because it has a 
meaningful Cramer-Rao bound and is more represen- 
tative of aircraft flight data. 

First we will look at the one-dimensional case 
where L$ is fixed at the correct value, because it 
is easier to grasp some of the characteristics of 
the cost function in one dimension. Then we will 
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look at the two-dimensional caae, where both Lp and 
Lj are varying. It ia important to remember that 
everything ahown in thia paper on coat functions ia 
baaed on computed tiee hiatories that are defined 
by Bq, (36). For every tiee hiatory we alight 
chooae (computed or flight data), a complete coat 
function ia defined. For the ca . of n variablea, 
the coat function definea a hyper-jurface of n ♦ 1 
dieenaiona. It night occur to ua that we could 
juat conatruct thia aurface and look for the mini- 
nua, avoiding the need to bother with the minimiza- 
tion algorithm, this ia not a reasonable approach 
because, in general, the number of variables is 
greater than two. Therefore, the cost function 
can be described mathematically but not pictured 
graphically. 

One-Dimensional Case 

To illustrate the many interesting aspects of 
cost functions, it is easiest to first look at cost 
functions having one variable. In an earlier sec- 
tion, the cost function of Lp and L^ was minimized. 
That cost function is most interesting in the Lp 
direction, therefore, the one-variable cost func- 
tion studied here is J(L p ). All discussions in 
this section are for J(Lp) with Lj equal to the 

correct value of 10. Figure 8 shows the cost func- 
tion plotted as a function of L p for the case where 

there is no measurement noise (no-noise case). As 
expected for this case, the minimum cost is zero 
and occurs at the correct value of Lp ■ -0.2500. 

It is apparent that the cost increases much more 
slowly for a more negative L p than fcr a positive 
L p . In fact, the slope of the curve tends to 
become less negative where Lp is more negative than 
-1.0. Physically this makes sense since the more 
negative values of L p represent cases of high 
damping, and the positive Lp represents an unstable 
system. Therefore, the p^ for positive L p becomes 
increasingly different from the measured time 
history for small positive increments in L p . For 
very large damping (very negative L p ). the system 
would show essentially no response. Therefore, 
large increases in damping result in relatively 
small changes in the value of J(L p ). 

In Fig. 9, the cost function based on the time 
history with measurement noise (noisy caae) is 
plotted as a function of L p . The correct value 
of Lp (-0.2500) and the value of L p (-0.3218) at 
the minimum of the cost (3.335) are both indicated 
on the figure. The general shape of the cost 
function in Fig. 9 is similar to that shown in 
Fig. 8. Figure 10 shows the comparison between 
the cost functions based on the noisy and no-noise 
cases. The comments relating to the cost function 
of the no-noise case also apply to the cost func- 
tion based on the noisy case. Figure 10 shows 
clearly that the two cost functions are shifted 
by the difference in the value of Lp at the mini- 
mum and increased by the difference in the minimum 
cost. One would expect only a small difference 
in the value of the cost when far from the mini- 
mum. This is because the "estimated* time history 
is so far from the measured time history that it 
becomes irrelevent as to whether the measured 


time history has noise added. Therefore, for large 
values of cost, the difference in the two cost func- 
tions should be small in comparison to the total 
cost. 

Figure 11 shows the gradient of J(Lp) plotted 
as a function of Lp for the noisy case. This is 
the function for which we were trying to find the 
zero (or equivalently, the minimum of the cost 
function) using the Gauss-Newton method that is 
discussed in a previous section. The gradient is 
zero at Lp - -0.3218, which corresponds to the 

value of the minimum of J(L p ). 

The difference between the Newton-Raphaon 
method (Eq. (14a)) and the Gauss-Newton method 
{Bq. 14b)) of minimization has been mentioned pre- 
viously. For this simple one-dimensional case, we 
can easily compute the second gradient both with 
the second term of Eq. (14a) (Newton-Raphson) , and 
without the second term (Gauss-Newton, Eq. (14b,). 
Figure 12 shows a comparison between the Newton- 
Raphson and the Gauss-Newton approximation second 
gradients. The Gauss-Newton second gradient 
(dashed line) always remains positive because it 
is the sum of quadratic terms (squared for the one- 
dimensional example). The Newton-Raphson second 
gradient can be positive or negative, depending 
upon the value of the second partial derivative 
with respect to L p . Other than the difference in 
sign for the more negative Lp, the two curves have 
similar shapes. 

As stated earlier, the O.'uss-Newton method can 
be shown to be superior to Newton-Raphson in cer- 
tain cases. We can demonstrate obvious cases of 
this with our example. An easy way to select a 
spot where problems with the Newton-Raphson method 
will occur is to look for places where the second 
gradient (slope of the gradient) is near zero or 
negative, Figure 11 has such a region near 
Lp = -1.0. If we choose a point where the gradient 

slope is exactly zero, we are forced to divide by 
zero in Eq. (12) with the Newton-Raphson method. 

This point is at Lp • -1.13 in Fig. 12, If the 
value of the slope of the gradient is negative, 
then the Newton-Raphson method will qo to very 
negative values of Lp. For very negative values 
of Lp, the cost becomes asymptotically constant and 
the gradient becomes nearly zero. In that region, 
the Newton-Raphson algorithm would diverge to nega- 
tive infinity. If aie slope of the gradient is 
positive but small, we still have a problem with 
the Newton-Raphson method. Figure 13 shows the 
first iteration starting from L p = -0.95 for both 
Gauss-Newton and Ne* ron-Raphson. The Newton- 
Raphson method selects a point where the tangent of 
the gradient at L p • -0.95 intersects the zero 
line. This results in the selection of an I. p of 
approximately 2.6 in the first iteration. From 
that value it requires many iterations to return to 
the actual minimum. On the other hand, the Gauss- 
Newton method selects a value for Lp of approxi- 
mately -0.09 and converges to the minimum to four- 
digit accuracy in two more iterations. With more 
complex examples a comparison of the convergence 
properties of the two algorithms becomes more 
difficult to visualize, but the problems are gener- 
alizations of the situation we have observed with 
the one-dimensional example. 



The usefulness of tha Cramer-Rao bound was dia- 
cusaad in tha Examples With Measurement Noise tac- 
tion- At this point it it uaaful to digrass 
briafly to discuss some of tha ramifications of tha 
Cramer-Rao bound for tha ons-dimens ions 1 casa- Tha 
Cramar-Rao bound only has manning for tha noisy 
case. In tha noiay example, tha astiaata of Lp ia 

-0.3218 and tha Cramer-Rao bound is 0.0579. Tha 
calculation of tha Cramar-Rao bound waa dafinad in 
tha pravioua aaction for both tha ona-diaansional 
and two-dimanaional axamplas. Tha Cramar-Rao bound 
is an astimata of tha standard daviation of tha 
astimata. Ona would axpact tha scattar in tha 
estimates of Lp to ^ of about the Sana magnitude 
as the estimate of tha standard daviation. For the 
one-dimensional case discussed hare, the range 
(Lp (-0.3218) plus or minus tha Cramar-Rao bound 
(0.0579) ) nearly includes the correct value of Lp 
(-0.2500). If noisy cases are generated for many 
time historias (adding different measurement noise 
to each time history), then the sample mean and 
sample standard daviation of the eatimatas for 
these cases can be calculated. Table 5 gives the 
sample mean, sample standard deviation, and the 
standard deviation of the sample mean (standard 
deviation divided by the square root of the number 
of cases) for 5, 10, and 20 cases. The sample 
mean, as expected, gets closer to the correct value 
of -0.2500 as the number of cases increases. This 
is also reflected by the decreasing values in 
column 3 of Table 5, which are estimates of the 
error in the sample mean. Column 2 of Table 5 
shows the sample standard deviations, which indi- 
cate the approximate accuracy of the individual 
estimates. This standard deviation, which stays 
more or less conttant, is approximately equal to 
the Cramer-Rao bound for the noisy case being 
studied here. In fact, the Cramer-Rao bounds for 
each of the 20 noisy cases used here (not- shown in 
the table) do not change much from the values found 
for the noisy case being studied. Both of these 
results are in good agreement with the theoretical 
characteristics lb of the Cramer-Rao bounds and 
maximum likelihood estimators in general. 

The examples shorn hare Indicate the value of 
obtaining more sample time histories (maneuvers). 
More samples improve confidence in the estimate of 
the unknowns. The same result holds true in ana- 
lyzing actual flight time histories (maneuvers)) 
thus it is always advisable to obtain several 
maneuvers at a given flight condition to improve 
the best estimate of each derivative. 

The size of the Cramer-Rao bounds and of the 
error between the correct value and the estimated 
value of Lp is determined to a large extent by the 

length of the time history and the amount of noise 
added to the correct time history. For the example 
being studied here, it is apparent from Fig. 5 that 
the asount of noise being added to the time history 
is large. The effect of the power of the measure- 
ment noise (GG*, Eqs- (3) and (4)) on the estimate 
of Lp for the time history is given in Table 6. 

The estimate of Lp is much improved by decreasing 

the measurement noise power. A reduction in the 
value of G to one-tenth of the value in the noisy 
example being studied yields an acceptable estimate 
of Lp. For flight data, the measurement noise is 

reduced by improving the accuracy of the output of 
the measurement sensors. 


Two-Dimensional Case 

In t. is section the cost function (dependant on 
both Lp and Lg) is studied. The no-noise case is 

examined first, followed by the noisy case. 

No-Noise Case . Even though the cost function 
is a function of only two unknowns, it becomes much 
more difficult to visualize than the one-unknown 
case. The cost function over a reasonable range of 
Lp and Lg is shown in Fig. 14. The cost increases 

vary rapidly in the region of positive Lp and large 
values of Lg. The reason is just an extension of 
tha argument for positive Lp given in the previous 
section. The shape of the surface can be depicted 
in greater detail if we examine only the values of 
tha cost function less than 200 for Lp less than 
1.0. Figure 15 shows a view of this restricted 
surface from the upper end of the surface. The 
minimum must lie in the curving valley that gets 
broader as we go to the far side of the surface. 

Now that we have a picture of the surface, we can 
look at the isoclines of constant cost on the 
Lp~versus-Lg plane. These isoclines are shown in 

Fig. 16. The minimum of the cost function is 
inside the closed isocline. The steepness of the 
cost function in the positive-Lp direction is once 
again apparent. Inside the closed isocline the 
shape is more nearly elliptical, indicating that 
the cost is nearly quadratic here, so fairly rapid 
convergence in this region would be expected. The 
Lp axis becomes an asymptote in cost as Lg 

approaches zero. The coBt is constant for Lg ” 0 
because no response would result from any aileron 
input. The estimated response is zero for all 
values of Lp, resulting in constant cost. 

Figure 16 shows the minimum value of the cost 
function, which, as seen in the earlier example 
(Table 2), occurs at the correct values for Lp and 
Lg of -0.2500 and 10, respectively. This is also 

evident by looking at the cost function surface 
shown in Fig. 17. The surface has its minimum at 
the correct value. As expected, the value of the 
cost function at the minimum is zero. 

Noisy Case . As shown before in the one- 
dimensional case, the primary difference between 
the cost functions for the no-noise and noisy 
cases was a shift in the cost function. In that 
instance, the noisy case was shifted so that the 
minimum was at a higher cost and a more negative 
value of Lp. In the two-dimensional case, the no- 
noise and noisy cost functions exhibit a similar 
shift. For two dimensions the shift is in both the 
Lp and Lg directions. The shift is small enough 

that the difference between the two cost functions 
is not visible at the scale shown in Fig. 14 or 
from the perspective of Fig. 15. Figure 18 shows 
the isoclines of constant cost for the noisy case. 
The figure looks much like the isoclines for the 
no-noise case shown in Fig. 16. The difference 
between Figs. 16 and 18 is a shift in Lp of about 
''.I. This is the difference in the value of Lp at 
the minimum for the no-noise and noisy cases. 
Heuristically, onu can see that the same would be 
true for cases with more than two unknowns. The 
primary difference^ between the two cost functions 
is near the minimum. 
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The next logical part of the coat function to 
examine ie near the minimum. Figure 19 ahova tha 
same view of the cost function for the noisy case 
as was shown in Fig. 17 for tha no-noise case. The 
shape is roughly tha same as that shown in Fig. 17. 
but the surface is shifted such that its minimum 
lies over Lp “ -0.3540 and Lg - 10.24, and is 

shifted upward to a cost function value of approxi- 
mately 3.3. 

To get a more precise idea of the coat of the 
noisy case near the minimum, we once again need to 
examine the iscolines. The isoclines (Fig. 20) in 
this regio.. are much more like ellipses than they 
are in Figs. 16 and 18. We can follow tha path of 
the minimization example used before by including 
the results from Table 4 on Fig. 20. The first 
iteration (L « 1) brought the values of Lp and L{ 
very close to the values at the minimum. The next 
iteration essentially selected the values at the 
minimum when viewed at this scale. One of the 
reasons the convergence is so rapid in this region 
is that the isoclines are nearly elliptical, demon- 
strating that the cost is very nearly quadratic in 
this region. If we had started the Gauss-Newton 
algorithm at a point where the isoclines are much 
leas elliptical (as in some of the border regions 
in Fig. 18), the convergence would have been much 
slower initially, but much the same as it entered 
the nearly quadratic region of the cost function. 

Before concluding our examination of the two- 
dimensional case, we need to examine the Cramer-Rao 
bound. Figure 21 shows the uncertainty ellipsoid, 

which is based on the Cramer-Rao bounds defined in 
an earlier section. The relationships between the 

Cramer-Rao bound and the uncertainty ellipsoid are 
discussed in Ref. 16. The uncertainty ellipsoid 
almost includes the correct value of Lp and Lj. 

The Cramer-Rao bound for Lp and Lj can be deter- 
mined from the projection of the uncertainty ellip- 
soid onto the Lp and L^ axes, and cong>ared with the 

values given earlier, which were 0.1593 and 1.116 
for Lp and L$, respectively. 

Estimation Using Flight Data 

In the previous several sections we examined 
the basic mechanics of obtaining maximum likelihood 
estimates from computed examples with one or two 
unknown parameters. Now that we have a grasp of 
these basics, we can explore the estimation of sta- 
bility and control derivatives from actual flight 
data. For the computationally much more difficult 
situation usually encountered using actual flight 
data, we will obtain the maximum likelihood esti- 
mates with the Illif f-Maine code (MMLE3 program) 
described in Pef. 17. The equations o* motion that 
are of interest are given in the Aircraft Equations 
of Motion section of this paper; the remainder of 
the equations are given in Ref. 17. 

In general, flight data estimation is fairly 
complex, and programs such as the Iliff-Maine code 
must usually be used to assist in the analysis. 
However, one must still be cautious about accepting 
the results; that is, the estimates must fit the 
phenomenology, and the match between the measured 
and computed time histories must be acceptable. 

This is true in all flight regimes, but one must be 
particularly careful in potential problem sit- 
uations such as ( 1 ) in separated flow at high na c h 


numbers or high angle of attack, (2) with unusual 
aircraft configurations , uch the oblique wing, 18 
or (3) with modern high-performance aircraft with 
high-gain feedback loops. In any of the above 
cases, one should be particularly careful where 
there are even small anomalies in the match. These 
anomalies may indicate ignored terms in the equa- 
tions of motion, separated flow, nonlinearities, 
sensor problems, or any of a long list of other 
problems • 

The following brief examples are intended to 
show how the above caveats and the computed e> . 
plea of previous sections can be used to ass : *■ ■. n 
the analysis. In the computed example, the Lr- 
ability of low-noise sensors, an adequate mociej , 
and several maneuvers at a given flight condition 
is shown. 

Hand Calculation Example 

Sometimes evaluation of a fairly complex flight 
maneuver can be augmented with a simple hand calcu- 
lation. One example of this can be found for the 
space shuttle. The space shuttle is a large 
double-delta-winged vehicle designed to enter the 
atmosphere from space and land horizontally. The 
entry control system consists of 12 vertical 
reaction-control-system (RCS> jets (six up-firing 
and six down-firing), 8 horizontal RCS jets (four 
left-firing and four right-firing), 4 elevon sur- 
faces, a body flap, and a split rudder surface. 

The locations of these devices are shown in Fig. 22. 
The vertical jets and the elevens are used for both 
pitch and roll control. The jets and elevons are 
used symmetrically for pitch control and asymmetri- 
cally for roll control. The space shuttle control 
system is described briefly in Ref. 6. 

The shuttle example used here is from a maneu- 
ver obtained at a Mach number of approximately 21 
and an angle of attack of approximately 40*. The 
controls being used for this lateral-directional 
maneuver are the differential elevons and the side- 
firing jets (yaw jets). The maneuver is shown in 
Fig. 23. Equations (15) to (31) describe the equa- 
tions of motion. A simplified approach can be used 
to determine some of the derivatives by hand. The 
approach is one that has been used since the begin- 
ning of dynamic analysis of flight maneuvers. In 
particular, for this maneuver e slope of the 
rates can be used to determine the yaw get control 
derivatives. This is possible for this example, 
even with a high-gain feedback system, because the 
yaw jets are essentially step functions and the 
slope of rates p and r can be determined before the 
vehicle and the differential elevon (aileron) 
responses become significant. The rolling moment 
due to yaw jet (Lyj) is particularly important for 
the shuttle 4n< j general, more difficult to 

obtain than the more dominant yawing moment due to 
yaw jet. Therefore, as an illustrative example, 
l YJ iB determined by hand. Figure 24 shows yaw jet 
and smoothed roll rate plotted at expanded scales. 
The equation for Lyj is given by 

Lyj “ pI x /(Number of yaw jets) (42) 

p - Ap/At - | ) ♦ ( 0 . 1 ) ( 43 ) 

Therefore, given that I x i 900,000 slug-ft 2 and 
the number of yaw jets is 4, Lyj s -2750 ft-lb. 




The same maneuver was analysed with KMLE3, and 
tha resulting natch ia shown in Fig. 25. Tha natch 
is vary good except for a snail nisnatch in p at 
about 6 sac. This snail nisnatch was studiad sep- 
arataly with MMLI3 and found to ba causad by a 
nonlinaarity in tha ailaron darivativa. Tha valua 
fron HHLI3 for Lyj is *2690 ft-lb, which for the 
accuracy usad hara ia assantially cha Sana valua at 
obtainad by tha sinplified mat hod. Tha ailaron 
darivativas would ba difficult to dataradna as 
accurataly as tha yaw jat darivativas. Although 
good astinatas can saldon ba obtainad with tha 
slopa nathod discussad hara, rough astinatas can 
usually ba obtainad to gain soma insight into 
valuas obtainad with HMLI3 (or any othar maxinun 
likalihood program). Thase rough astinatas can 
than ba usad to halp ai.plain unaxpactad valuas of 
astinatas fron an astination program. 


convarganca ara dasirabla, but not sufficiant, 
conditions for a succassful analysis. Flgura 27(a), 
although not parfact, is a vary good natch. Tha 
convarganca can bast ba avaluatad by looking at tha 
nornalisad cost in tha last row of Tabla 7. Tha 
cost has rapidly and nonotonically convargad in 
four ltarations, and it ranains at tha convargad 
cost. Thasa factors ara convincing avidanca that 
tha convarganca in conplata. Tharafora, tha crl- 
taria of natch and convarganca ara satisfiad in our 
axanpla. In soma casas wa night ancountar cost 
that doas not convarga rapidly (in four to six 
ltarations) or nonotonically, or stay "exactly" at 
tha nininun valua. Thasa situations usually indi- 
cata at laast a snail problem in tha analysis. 

Thasa problems, if found, are usually traced to 
a data problem, an inadequate mathematical modal, 
or a maneuver that contains a marginal amount of 
infornation. 


Sometimes a flight example becomes too complex 
to get anything other than qualitative estimates by 
hand. Tha determination of the rudder derivative 
for tha F-8 aircraft with the yaw augmentation sys- 
tem on demonstrates this difficulty. Figure 26 
shows an example of this difficulty for tha F-8 
data. This example, taken from Ref. 19, includes 
an aileron pulse and a rudder pulse. Although an 
independent pilot .udder pulse is input during tha 
maneuver, tha rudder is largely responding to the 
lateral acceleration feedback. When the rudder is 
moving, several other variables are also moving, 
thus making it difficult, to use the simplified 
approach just discussad. However, C n j can be 

roughly determined when the rudder moves, approxi- 
mately 1.7 sec from the start of the maneuver. 

Most of the slope of yaw rate is caused by the 
rudder, but a poor estimate would be obtained using 
the hand calculation. 

Cost Function for Full Aircraft Problem 

The analysis of a lateral-directional maneu- 
ver obtained in flight typically has from 15 to 25 
unknown parameters (as shown in Eqs. ( 15 ) and 
(31)), in contrast to the one or two in the simple 
aircraft example- This makes detailed examples 
unwieldy and any graphic presentation of the cost 
function impossible . Therefore, in this section we 
are primarily examining the estimation procedure 
and the process of the minimization. 

For our flight example, we have chosen a 
lateral-directional maneuver, with both aileron and 
rudder inputs, that has 17 unknown parameters. The 
data are from the oblique wing aircraft 18 with the 
wing unskewed during the maneuver. This example 
was chosen because it is a typical maneuver. The 
tine history of the data and the subsequent output 
of MMLE3 have been published in Ref. 20. The tabu- 
lar results of the analysis are shown in Table 7. 
The natch between the measured tine history (solid 
lines) and the estimated (calculated) time history 
(dashed lines) is shown as a function of iteration 
in Fig. 27. Figures 27(a) to (e) are for itera- 
tions 0 to 4, respectively. Table 7 shows that the 
cost remains unchanged after four iterations. A 
similar result was obtained for the two-dimensional 
simple aircraft example in Fig. 6 and Tabla 4. 

Of the many things the analyst must consider 
in obtaining estimates, the two most important 
ones are how good is the match and how good Is the 
convergence. A satisfactory match and monotonic 


Table 7 also shows that the startup values of 
all the coefficients are zero for the cctrol and 
bias variables. Hind tunnel estimates could have 
been used for starting values, but the convergence 
of the algorithm is not very dependant on the 
startup values. As part of the startup algorithm, 
the MMLE3 program normally holds the derivatives of 
the state variables constant until after the first 
iteration, as is evident in Table 7. 


Figure 27(a) shows the match between the 
measured and computed data for the startup values. 
The match ia very poor because the startup values 
for the control derivatives are all zero, so the 
only motion is in response to the initial condi- 
tions. The control derivatives and biases are 
determined on the first iteration, resulting in 
the much improved match shown in Fig. 27(b). The 
match after two iterations, shown in Fig. 27(c), 
is improved as the program further modifies the 
control derivatives and, for the first time, 
adjusts the derivatives affecting the natural 

frequency (C-„ and Co_). By the third iteration 
"B *B 

(Fig. 27(d)), th .provement in the match is 

almost complete, because minor adjustments to the 
frequency are made and the damping derivatives are 
changed. Fig. 27(e) shows the match when all but 
the most minor derivatives have ceased to change. 


Several general observations can be made based 
on this well behaved example. The strong or most 
important coefficients have essentially converged 
in three iterations. The same effect was seen in 
the simple example - that is, Lj converged faster 


than Lp (Table 4). Some of the less important or 

second-order coefficients have only converged to 
two places after three iterations and are st ' 11 
changing by one digit in the fourth place at the 
end of six iterations- Another observation is that 

for some coefficients (C, , C n , , and C«, ) even 

r “4 a o r 

though the sign is wrong after the first iteration, 
the algorithm quickly selects their correct values 
once the important derivatives have stabilized. 


In general, if the analysis of a maneuver has 
gone well, we do not need to spend much time inspec- 
ting a table analogous to Table 7. However, if 
there have been problems in convergence or in the 
quality of the fit, a detailed inspection of such a 
table may be necessary. The data may show an impor- 
tant coefficient going unstable at an early itera- 
tion, which could cause problt a later. If the 
starting values are grossly in error, the algorithm 
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is driven a long way from reasonable values and 
than for many reasons does not behave well • Occa- 
sionally the algorithm alternately selects from two 
diverse sets of values of two or more coefficients 
on successive iterations, behaving as if the shape 
of the cost function were a narrow multidimensional 
valley analogous to but more extreme than the two- 
dimensional valley shown in Pigs, 18 and 20, 

Cramer-Rao Bounds 


The earlier sections regarding the computed 
example have shown that the Cramer-Rao bound is 
a good indicator of the accuracy of an estimated 
parameter. The Cramer-Rao bounds can be used in 
a similar, but somewhat more qualitative, fashion 
on flight data. The Cramer-Rao bounds that are 
Included in MMLE3 (as well as many other maximum 
likelihood estimation programs) have been useful 
in determining whether estimates are good or bad. 
The aircraft example discussed here has been 
reported previously (for example, in Refs, 1 and 
16), However, this example of the use of the 
Cramer-Rao bound in the assessment of flight- 
derived estimates is pertinent to the thrust of 
this paper. Figure 28 shows estimates of C n as 


a function of angle of attack for the PA-30 twin- 
engine general aviation aircraft 21 at three flap 
settings. There is a significant amount of scat- 
ter, which makes the reliability of the information 
on C n questionable. The data shown are the esti- 


mates from the MMLE3 program, which also provides 


the Cramer-Rao bounds for each estimate. Past 


1 A 

experience* has shown that if the Cramer-Rao bound 
is multiplied by a scale factor (the result some- 
times called the uncertainty level 1,16 ) and plot- 
ted as a vertical bar with the associated esti- 
mate, it helps in the interpretation of flight- 
determined results. Figure 29 shows the same data 
as Fig. 28, with the uncertainty levels now inclu- 
ded as vertical bars. The estimates with small 


uncertainty levels (Cramer-Rao bounds) are the best 
estimates, as was discussed earlier in the section 
on Cramer-Rao bounds for the one-dimensional case. 
The fairing shown in Fig. 29 goes through the esti- 
mates with small Cramer-Rao bounds ar.' ignores the 
estimates wit ~.arge bounds. One can have great 
confidence in the fairing of the estimates, because 
the fairing is well defined and consistent when the 
Cramer-Rao bound information is included. In this 
particular instance, the estimates with small 
bounds were from maneuvers where the aileron forced 
the motion, and the large bounds were from maneu- 
vers where the rudder forced the motion. There- 
fore, in addition to aiding in the fairing of the 

estimates, the Cramer-Rao bounds help show that the 
aileron-forced maneuvers are superior for esti- 
mating C n for the PA-30 aircraft. 


This example illustrates that the Cramer-Rao 
bounds are a useful tool in assessing flight- 
determined estimates, just as they were found use- 
ful for the simple aircraft example with computed 
data. 


example showed the advantage of low measurement 
noise, multiple estimates at a given condition, and 
the Cramer-Rao bounds, and the quality of the match 
between the measured and computed data. The flight 
data showed that many of these concepts still hold 
true even though the dimensionality of the cost 
function makes it impossible to plot or visualise. 
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Table 1 Values of computed time 
history with no measurement ncise 


i 

6, deg 

p, deg/sec 

1 

0 

0 

2 

1 

0.9754115099857 

3 

1 

2.878663149266 

4 

1 

4.689092110779 

5 

1 

6.411225409939 

6 

1 

8.049369277012 

7 

1 

9.607619924937 

8 

0 

10.11446228200 

9 

0 

9.621174135646 

10 

0 

9.151943936071 


Table 2 Pertinent values as a function of iteration 


L 

Lp(L) 

Ci(L) 

♦ CL) 

4>(L) 

•*L 


0 

-0.5000 

15.00 

0.9048 

2.855 

21.21 


1 

-0.3005 

9.888 

0.9417 

1 .919 

0.5191 


2 

-0.7475 

9.996 

0.9517 

1.951 

5.083 x 

10-4 

3 

-0.2500 

10.00 

0.9512 

1.951 

1.540 x 

10- 9 

4 

-0.2500 

10.00 

0.9512 

1.951 

1.060 x 

10- 1 * 


Table 3 
tory with 

Values 

added 

of computed time his- 
measurement noise 

Table 

4 Pertinent values 

as a function of 

iteration 

i 

{, deg p, deg/sec 

L 

Lp(L) 

LfiCL) 

♦ CL) 

*( L) 

Jt 

1 

2 

n 

n 



u 

1 

0.4875521781881 

0 

-0.5000 

15.00 

0.9048 

2.855 

30.22 

3 

1 

3.238763570696 

1 

-0.3842 

10.16 

0.9260 

1.956 

3.497 

4 

1 

3.429117357944 

2 

-0.3518 

10.23 

0.9321 

1 .976 

3.316 

0 

1 

6.286297353361 

3 

-0.3543 

10.25 

0.9316 

1 .978 

3.316 

6 

1 

6.953798550097 

4 

-0.3542 

10.24 

0.9316 

1 .978 

3.316 

7 


10.80572930119 

5 

-0.3542 

1C. 24 

0.9316 

1.978 

3.316 

8 

0 

9.739367269447 








9 

10 


9.788844525490 

7.382568353168 


Table 5 Mean and standard deviations for estimates of Lp 


Number of 
cases, N 

sample mean, 

u<V 

Sample standard 
deviation, a(Lp) 

Sample standard 
deviation of the 

mean, alLp)//W~ 

5 

-0.2668 

0.0739 

0.0336 

10 

-0.251 . 

0.0620 

0.0196 

20 

-0.2452 

0.0578 

0.0129 


12 


\ — 


4J?*: 



Table 6 Eatieat. of Lp and Crae.r-Rao bound aa 
a function of the aquare root of noiaa power 


Square root of 
noiaa power, G 

Batleate 

of Lp 

Craeer-Rao 

bound 

0.0 

-0.2500 


0.01 

-0.2507 

0.00054 

0.05 

-0.2535 

0.00271 

0.10 

-0.2570 

0.00543 

0.2 

-0.2641 

0.0109 

0.4 

-0.2783 

0.0220 

0.8 

-0.30T1 

0.0457 

1 .0 

-0.3218 

0.0579 

2.0 

-0.3975 

0.1248 

5.0 

-0.6519 

0.3980 

10.0 

-1.195 

1.279 


Table 7 Stability and control derivatives as a function of iteration for flight maneuver 


Iteration 




Iteration 




L 

1 

2 

3 

4 

5 

6 

C *0 

- 0.008500 

- 0.008500 

- 0.007959 

- 0.008347 

- 0.008375 

- 0.008364 

- 0.008364 

c *e 

- 0.0002500 

- 0.0002500 

- 0.0003141 

- 0.0003580 

- 0.0003572 

- 0.0003571 

- 0.0003571 

Cn e 

0.001000 

0.001000 

0.001159 

0.001243 

0.001230 

0.001230 

0.001230 

% 

- 0.2500 

- 0.2500 

- 0.3393 

- 0.3584 

- 0.3581 

- 0.3581 

- 0.3581 

Cn P 

- 0.02500 

- 0.02500 

- 0.04356 

- 0.04537 

- 0.04512 

- 0.04599 

- 0.04600 

C *r 

- 0.05000 

- 0.05000 

0.06790 

0.07044 

0.06972 

0.06973 

0.06974 

C "r 

- 0.00800 

- 0.08000 

- 0.1327 

- 0.1033 

- 0 . 106 b 

- 0.1062 

- 0.1062 

C ‘*a 

0 

0.0008009 

0.001000 

0.001067 

0.001069 

0.001069 

0.001069 

Cn «a 

0 

- 0.00004604 

6.786 » 10" 7 

0.00001129 

0.00001096 

0.00001068 

0.00001069 

c v« t 

0 

0.005935 

0 . 0 C 2064 

0.001456 

0.001546 

0.001548 

0.001548 

C ^r 

0 

- 0.00005068 

0.00005764 

0.0001043 

0.0001059 

0.0001055 

0.0001055 

Cn »r 

0 

- 0.0007329 

- 0.0009333 

- 0.0008875 

- 0.0008972 

- 0.0008961 

- 0.0008961 

c v 0 

0 

- 0.05109 

- 0.02691 

- 0.02362 

- 0.02420 

- 0.02419 

- 0.02420 

C Y 0 + ®0 

0 

- 0*03370 

- 0.01370 

- 0.01 117 

- 0.0115 

- 0.01156 

- 0.01156 

c *o 

0 

- 0.0007096 

- 0.001629 

- 0.002021 

- 0.002031 

- 0.002028 

- 0.002028 

c "o 

0 

0.005864 

0.007300 

0.007140 

0.007175 

0.007169 

0.007169 

*0 

0 

0.2121 

0.1626 

0.1482 

0.1506 

0.1506 

0.1506 

J/(N - 1) 

731 .5 

65.00 

11.23 

4.8265 

4.701 

4.701 

4.701 
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Damping In roll, L p 


Pig. 9 Coat function as a function of Lp for noisy 

casa. 



Pig. 11 Gradient of J(L p ) as a function of Lp 
for. noisy case. 
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Damping in roll, L p 

Pig. 10 Comparison of the coat functions for 
the no-noise and noisy cases. 



Pig. 12 Comparison of Newton-Raphson and Gauss- 
Newton values of the second gradient for the noisy 
case. 



Fig. 13 Comparison of first iteration step size 
for the Nevton-Raphson and Gauss-Newton algorithms 
for the noisy case. 
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Fig. 24 Bitemples of obtaining Lyj by simple calcu- 
lations for the shuttle data from Fig . 23. 
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Fig, 23 Lateral-directional space shuttle maneuver 
at a Mach number of 21. 
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Fig. 25 MMLE3 match of maneuver shown in 
Fig. 23. 
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Fig. 26 Lateral-directional maneuver 
from F-8 aircraft with augmentation on. 
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(a) Zero iterat ion. 

Fig . 27 Match between measured and computed 
time histories as a function of iterat ion . 
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(b) One iteration. 
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Fig. 27 Continued . 
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(d) Three iterations . 
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(e) Four iterations . 


Fig. 27 Concluded. 



